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Abstract 


A numerical method for the quasi-neutral two-fluid (QNTF) plasma model is described. The basic equations are ion 
and electron fluid equations and the Maxwell equations without displacement current. The neglect of displacement 
current is consistent with the assumption of charge neutrality. Therefore, Langmuir waves and electromagnetic waves 
are eliminated from the system, which is in clear contrast to the fully electromagnetic two-fluid model. It thus reduces 
to the ideal magnetohydrodynamic (MHD) equations in the long wavelength limit, but the two-fluid effect appearing 
at ion and electron inertial scales is fully taken into account. It is shown that the basic equations may be rewritten in a 
form that has formally the same structure as the MHD equations. The total mass, momentum, and energy are all written 
in the conservative form. A new three-dimensional numerical simulation code has been developed for the QNTF 
equations. The HLL (Harten-Lax-van Leer) approximate Riemann solver combined with the upwind constrained 
transport (UCT) scheme is applied. The method was originally developed for MHD (Londrillo & Del Zanna, 2004!), 
but works quite well for the present model as well. The simulation code is able to capture sharp multidimensional 
discontinuities as well as dispersive waves arising from the two-fluid effect at small scales without producing V • B 
errors. It is well known that conventional Hall-MHD codes often suffer a numerical stability issue associated with 
short wavelength whistler waves. On the other hand, since finite electron inertia introduces an upper bound to the 
phase speed of whistler waves in the present model, our code is free from the issue even without explicit dissipation 
terms or implicit time integration. Numerical experiments have confirmed that there is no need to resolve characteristic 
time scales such as plasma frequency or cyclotron frequency for numerical stability. Consequently, the QNTF model 
offers a better alternative to the Hall-MHD or fully electromagnetic two-fluid models. 

Keywords: collisionless plasma, magnetohydrodynamics. Hall magnetohydrodynamics, HLL Riemann solver, 
constrained transport 


1. Introduction 

Understanding of a rich variety of nonlinear phenomena in space, astrophysical, and laboratory plasmas requires 
numerical simulations at various levels of approximations. At the largest scale, magnetohydrodynamic (MHD) de¬ 
scription is useful because the scale-free nature of MHD allows us to conduct simulations with a realistic scale size. 
On the other hand, physics at kinetic scales (i.e., ion and electron inertial lengths) must be taken into account in cases 
where it plays the key role. A well-known example is the diffusion region of collisionless magnetic reconnection, 
in which the kinetic effect is essential for violating the frozen-in condition. Fully kinetic particle-in-cell (PIC) sim¬ 
ulations have been used to investigate such problems. It is important to point out that the characteristics of spatially 
localized tiny regions may ultimately affect even the global dynamics of the system. 

Although it is believed that physics beyond MHD ultimately needs to be incorporated properly even for the mod¬ 
eling of macroscopic phenomena, it is still a formidable task, albeit not impossible, to perform fully kinetic PIC 
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simulations at a macroscopic scale. In practice, it is desirable to start with a simpler model and gradually proceed to¬ 
ward better (but more complicated) physics models with less approximations. Hall-MHD is one of such better models 
in the sense that it takes into account physics at the ion inertial scale. The hybrid simulation model that deals with 
kinetic ions and a massless fluid electron can be considered as a “kinetic version” of Hall-MHD. The Hall-MHD and 
hybrid models are therefore believed to be possible alternatives to MHD for the next generation global modeling. 

Although both Hall-MHD and hybrid have been well established standard models for simulations of collisionless 
plasmas, it is well known that they often suffer a numerical difficulty due to high frequency whistler waves. The dis¬ 
persion relation of whistler waves to oc k 2 leads to the increase in the phase speed at short wavelength without bound. 
This is generally thought of as a source of numerical instability. A common strategy to stab ilize such simulations is to 
introd uce ad-hoc numerical dissipation such as hyper-resistivity in the code ( Ma & Bh attacharieS. 200 it Shav et all 
2001). However, it is not easy to control the amount of numerical dissipation with this kind of approach. Furthermore, 
since the strategy is quite different from the philosophy of modern high-order shock capturing schemes, this makes it 
difficult to extend such codes to the Hall-MHD regime. Although one may use an impli cit sc heme to circumvent the 


problem, this will make implementation of the algorithm much more complex (e.g., Arnold et al 

!2008l). 


2008; Toth et al 


A more straightforward approach is to employ the fully electromagnetic two-fluid (EMTF) plasma model in 
which the full set of Maxwell equations are coupled with two separate (i.e., ion and electron) fluids equations 


( Shumlak & Loverichl 2003 : Loverich & Shumlak , 2005 : Hakim et al. . 2006 : Srinivasan & Shumlaki 201 it Kumar & Mishral 


20121) . The phase speed of whistler waves has an upper bound in this system due to the presence of finite electron iner¬ 
tia. On the other hand, since it is essentially a fluid counterpart of the PIC simulation, it must deal with high frequency 
Langmuir waves as well as electromagnetic waves. Numerical stability requires that these waves should ade quat ely 
be resolved by the simulation time step unless more complicated implicit schemes are employed (IKumar & Mishral 
2012). In general, however, these high frequency waves are not of interest as far as macroscopic dynamics is con¬ 
cerned. The neglect of displacement current (which implies the charge neutrality) in the Maxwell equations is indeed 
a reasonable assumption if one considers non-relativistic problems, although this does not in general apply to highly 
relativistic plasmas (e.g., Amano & Kirk, 2013). 

There is also a way to incorporate finite electron inertia effect into Hall-MHD/hybrid without resorting to the 
full set of Maxwell eq uations. Conventionally , finite electron inertia effect has been included as a correction to 
the magnetic field (e.g., Kuznetsova et all 1998 : Shav et all 1998 . 2001 : Nakamura et all 2008). Although most of 
previous studies adopted some kind of simplification, the finite electron inertia effect if appropriately included can 
correctly introduce an upper bound to the phase speed of whistler waves. On the other hand, the motivation for 
these studies was to initiate spontaneous magnetic reconnection without relying on an anomalous resistivity model. 
Therefore, a possible advantage of finite electron inertia effect on the numerical stability issue has not been paid much 
attention. We have recently shown that, by modifying the procedure to incorporate finite electron inertia into the 
model, hybrid simulations can be made much more robust particularly in low-density regions where whistlers become 
problematic ( Amano et all 2014 1. This was made possible by implementing finite electron inertia as a correction to 
the electric field (i.e., the generalized Ohm’s law), which is then used to update the magnetic field. This physically 
more consistent approach gives a natural way to handle even a pure vacuum region in a hybrid code. It is quite natural 
to expect that essentially the same methodology can be applied to Hall-MHD equations, because kinetic ion dynamics 
does not play a role for dispersion of whistler waves. 

In the present paper, we consider a system consisting of two-fluid equations coupled with Maxwell equation 
without displacement current (Darwin approximation), which we call the quasi-neutral two-fluid (QNTF) model. As 
we will see in the next section, it is approximately the same as the Hall-MHD equations with finite electron inertia, but 
terms dropped in previous studies are retained for consistency. Consequently, the total mass, momentum, and energy 
including finite electron contributions are all written in the conservative form. The conservation laws coupled with the 
induction equation for the magnetic field have the same formal structure as the MHD equations, which thus may be 
solved by a known conservative scheme. Because of the neglect of displacement current, there are no high frequency 
waves such as Langmuir or electromagnetic waves, and the number of eigenmodes is indeed the same as MHD. The 
system correctly reduces to the ideal MHD in the long wavelength limit. Therefore, we think that it provides a natural 
extension of MHD having desirable properties both in terms of physics and numerics. 

We have developed a three-dimensional (3D) numerical simulation code to solve the proposed system of equations. 
We employ the single-state HLL (Harten-Lax-van Leer) approximate Riemann solver as a building block. It only 
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requires the maximum characteristic speed, and is independent of detailed information on the eigenmode structure. 
The scheme is thus suitable to the QNTF equations because eigenmode decomposition for this system should certainly 
be much more laborious task than for MHD. In addition, we adopt the Upwind Constrained Transport (UCT) scheme 
to keep the divergence error of the magnetic field within mac hine accurac y (Londrillo & Del ZannaL 120041) . The UCT 
scheme is based on the Constrained Transport (CT) scheme ( Evans & Hawlevl 19881) . but is designed specifically to 
be consistent with an underlying Riemann solver. Although it was originally developed for MHD, we found it is 
useful for the QNTF equations as well. With these numerical techniques, our simulation code is able to capture sharp 
discontinuities as well as dispersive waves arising from the two-fluid effect at the same time even in multidimensions 
without violating the divergence-free property. 

In the next section, we introduce the QNTF model. The characteristics of the model and its advantages over the 
Hall-MHD and EMTF models are discussed. Section 3 is devoted to numerical algorithm used in our code. Numerical 
results of several benchmark problems are discussed in section 4. Finally, conclusions are given in section 5. 


2. Quasi-neutral Two-fluid Model 

2.1. Basic Equations 

We start with the following fluid equations for a particle species s (i and e for ion and electron respectively) of 
charge q s and mass m s : 


-p s + W-(p s y s ) = 0 


4 PsVs + V ■ (PvVvV, + Ps I) = — Ps (e + — X b) 
ot m s \ c / 


Ps V 


L P. s v ■ E, 


( 1 ) 

( 2 ) 

(3) 


where p s , v s , p s are the mass density, bulk velocity, and (scalar) pressure (with I being the unit tensor), respectively. 
We here assume a polytropic equation of state with a specific heat ratio denoted by y (independent of particle species). 
The right-hand side of the above equations represents the Lorentz force with the electromagnetic field E, B, and the 
speed of light c. 

Since we only consider low frequency phenomena, the displacement current in the Maxwell equations is ignored. 


B = -VxE 

c dt 

An 

V x B = — J. 

c 


(4) 

(5) 


As usual, the electric current density is given by a sum of contributions from ions and electrons 


T q i qe 

J — Pi^i Pe^e 

mi m (> 


( 6 ) 


and we always assume charge neutrality 


0 = 


qt 



(7) 


which may also be written as n, = n e for q, = -q e = e. Here, n s is the number density and e is the elementary 
charge. This is indeed consistent with the neglect of the displacement current, because the longitudinal part of the 
displacement current represents charge-density fluctuations. The solenoidal condition for the magnetic field 


VB = 0 


( 8 ) 


gives a constraint that must be satisfied. 
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The above equations have source terms in the right-hand side and the energy and momentum are not strictly 
conserved quantities in numerical simulations. Instead, we use the following computationally more convenient con¬ 
servative form of equations that are obtained by taking sum of the two species: 

d 

—U + V ■ F = 0, (9) 

at 


where 


U = 


V 2 


PP>1 


Pi + Pe 

ppfi + P<Ae 

1 2 1 
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( 10 ) 


and 
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p,V,V; +p e \ e Ve + Pi + Pe + 
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( 11 ) 


represent conservative variables and their corresponding fluxes. Here Eqs. (j4]|7} have been used to rewrite the Lorentz 
force on the right-hand side into the above conservative form. Note that the same strategy was recently used in fully 
relativistic two-fluid simulations^.e., relativistic version of EMTF), so that the total energy and momentum become 
strictly conserved quantities (Amano & Kirk-, 120131) . 

As shown in Appendix A the generalized Ohm’s law for the present system may be written as 


(A + c 2 V x Vx) E = -- x B + V ■ n + 77 A J. 


( 12 ) 


Here we have introduced resistivity 77 in a rather ad-hoc manner to take into account phenomenological collisions, 
although it is absent (77 = 0) in the original ideal two-fluid equations. The moment quantities A, T, II appearing in the 
above equation are defined as follows 


A a Qe 2,2 

= 4 t: pi— + An Pe — = oj ■ + oj 

m m- ' F 

l c 

(13) 

q 2 q 2 e 7 7 

= Anpi—Mi + Anp e —\ e = 0J pl \, + to pe \ e 

tfl ■ 

l C 

(14) 

Anq, A nq e 

= - (P/V/V/ + pj) + - (pe'/e'/e + Pe I) , 

(15) 


where io 2 s = Anp s q 2 /m 2 is the plasma frequency for a particle species s. The connection between Eq. ( 1 1 2b and 
conventional Ohm’s laws will be discussed in the next subsection. We should emphasize that, aside from the resistivity 
77 , this form of the generalized Ohm’s law is exact. It is obtained from the basic equations without any approximations 
or assumptions. Notice that the above equation is an implicit equation for the electric field. Therefore, in general, 
matrix inversion is needed to obtain the electric field. Detail will be discussed later in section [3~4l 

Eqs. 03 clearly indicate that the density and velocity for each species are not independent. Although the electron 
density and velocity appear frequently in this paper for clarity of notation, they must always be replaced by 


qi/rrii 

Pe ~ Pi . 

q e /m e 

111! C 

\ e = v/--—V x B 

qt Anp { 


(16) 

(17) 
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in actual calculations. 

In addition, the ion to electron temperature ratio denoted by r = T, / T e is assumed to be a given constant throughout 
in this paper. This implies that the energy exchange between the species occurs instantaneously and the temperature 
ratio quickly relaxes to the prescribed constant value r (typically chosen to be t = 1). We adopt this simplification 
because the energy distribution among different species through a dissipative process (such as a shock wave) is not 
known a priori within the framework of a fluid model. It crucially depends on complicated kinetic physics in an 
unresolved dissipative layer. We note that there is no fundamental difficulty in dealing with T, and T e independently 
without such a simplifying assumption. Indeed, in such a case, we have found that there appears an entropy-like mode 
which is a pressure-balanced structure across which only ion and electron temperatures are exchanged with keeping 
the total gas pressure (and other quantities) unchanged. However, in general, it is only the total gas pressure behind 
the dissipative layer that we can correctly predict (from the Rankine-Hugoniot relations), and the temperature ratio is 
likely to be affected by numerical dissipation. Since the total gas pressure does not depend on a particular choice of r, 
we have not observed any noticeable differences between the simplified and the more rigorous implementations. We 
thus think the assumption adopted here is a reasonable simplification. 

In summary, the fluid quantities and magnetic field are updated respectively by using Eq. (0 and Eq. 0. The 
electric field appearing in these equations is determined by the generalized Ohm’s law Eq. (IT2l) . These equations and 
the relationship Eq. 00 with the constant ion and electron temperature ratio r close the system of equations, which 
is called the QNTF model. It is important to mention that the number of eigenmodes in this system is seven, which is 
the same as MHD. Namely, one may consider u = {/>,, v,-, p,, B} as primitive variables (with the V • B = 0 constraint). 
Given the ion quantities and the magnetic field, the electron quantities are automatically determined by Eq. (1161117b . 
As in the case of MHD, the electric field is also essentially a dependent variable, since it is completely specified by 
the primitive variables. 

Obviously, the advantage of using the conservative form Eq. 0 instead of the original two-fluid equations written 
separately is that the exact conservation of total energy and momentum may be guaranteed if a conservative scheme 
is used for numerical computation. Furthermore, since the QNTF equations in the conservative form are very similar 
to MHD equations, one may use numerical methods developed for MHD with relatively minor modifications. In 
particular, the absence of the Lorentz force as the source term gives an advantage, because otherwise it would possibly 
impose a constraint on the time step for numerical stability. As we will demonstrate with the numerical examples 
presented in section0 the time step of our code is restricted by the fastest wave mode, and there is no need to resolve 
characteristic time scales such as the cyclotron frequency. 

For the sake of completeness, we here give explicit formulae to calculate the primitive variables from the conser¬ 
vative variables. Given the conservative variables U = {!). M, K), one can calculate p,- and v,- according to 



( 18 ) 


(19) 


The electron density and velocity are determined by using Eqs. (1161117b . from which one may obtain the ion pressure 
Pi (and also electron pressure p e = pJt) as follows 



( 20 ) 


2.2. Model Characteristics 

It is helpful to introduce approximations to the generalized Ohm’s law Eq. (IT2l) in a step-by-step manner to see 
the relationship with more familiar forms of Ohm’s law. First, notice that A and T are the sum of contributions of the 
two species that are proportional to the plasma frequency. Similarly, the contributions to II are inversely proportional 
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to the mass. Therefore, one may safely ignore the ion fluid contributions and adopt the following approximation: 


A ~ <4e> 
r ~ u>p e Y e , 

r, 4nq e 

n » - (Pe'fe'/e + Pel) , 

m e 


( 21 ) 

( 22 ) 

(23) 


for nig/m, «: 1. This yields the following Ohm’s law: 

2 \ 

C V 171 

1 + -=-V X VX E = -— X B - — V • (p e V e V e + p e I) + pi, (24) 

Upe J C Pe 

which includes correction terms approximately describe the finite electron inertia effect. Finite electron inertia codes 
that have been used in previous studies, adopt essentially the same approximation (except for some minor differences). 
Now, it is easy to understand that the second term in the left-hand side of the above equation is smaller than the first 
term by a factor k 2 c 2 /a) 2 pe , which thus can be ignored for scale length longer than the electron inertial length c/a> pe . If 
we further drop the convective derivative term (oc p e \ e y e ), we obtain Ohm’s law for the Hall-MHD regime (i.e., with 
a massless electron fluid): 


E = —-xB + — — (VxB)xB--V/7 e + //J, (25) 

c q t 4npi p e 

where the electron velocity has been replaced by Ea. (fI7l i. It is well known that we obtain Ohm’s law for the MHD 
regime by taking the long wavelength (longer than the ion inertial length), and cold electron limit ( p e —♦ 0). This 
analysis confirms that Eq. (IT2l > indeed generalizes the known forms of Ohm’s law. From the rigorous Ohm’s law 
without approximations, we see that the magnetic field convection (or frozen-in) velocity must be given by 


r 

A 


(^pe y e 


■ u 2 pi Vi 


OJZ 


■ ^>n 


+ — (V/ - Ye) + O 
171; 


\ 


(26) 


rather than the electron velocity v e . The finite electron inertia correction appearing in the above equation has been 
ignored in previous studies, although one may expect it to be a relatively minor correction. 

There is another concern associated with the approximation Eq. (l24l) . Namely, it will break Galilean invariance 
due to the appearance of an unphysical electric field. To see this, consider for simplicity a cold (p, = p e = 0), current- 
free (V x B = 0) MHD flow. From Eq. ifTTb . it follows that the ion and electron flow velocities must be the same 
Vi = \ e = V. In this case, the ion and electron contributions to II in the exact equation cancel with each other: 


n = 


4nqi 
- Pi + 


nij 



VV = 0 


(27) 


because of the charge neutrality assumption Eq. 0. On the other hand, it is obvious that the approximate expression 
of II remains finite and produces an unphysical electric field when the flow speed or density has spatial variation (i.e., 
V • II + 0). Although the magnitude of the electric field will be small in comparisons with other contributions unless 
the flow speed is highly supersonic, it is better to keep the ion contributions included for consistency. The problem 
arises because careless neglect of terms of order 0(m e /mj ) breaks the local momentum conservation law. 

In the present model, the finite electron inertia effect is fully taken into account in the sense that it is correct to 
all orders of This enables us to write the equations for total energy and momentum including the electron 

contributions in the conservative form. This property has been missing in existing finite electron inertia codes. Con¬ 
sequently, the model is valid even for a pair plasma m, = m e , although application of a non-relativistic pair plasma 
model would be limited in practice. This may be sometimes useful for a control experiment because in this case the 
Hall term disappears and the perfect symmetry is preserved as in the case of MHD. 

One may recognize the present model as a better alternative to the Hall-MHD (with or without finite electron 
inertia) or EMTF models. In Hall-MHD, dispersive whistler waves often pose numerical difficulty because the phase 
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speed increases without bound at short wavelength. By taking into account finite electron inertia effect, there appears 


an upper bound in the phase speed that may improve numerical stability (see Appendix B for detail). The fact that the 
basic equations are written in the conservative form makes it easy to apply a known scheme for a hyperbolic conser¬ 
vation law. As we will see in the numerical examples discussed in section [4] the QNTF model automatically reduces 
to MHD in the long wavelength limit. This property allows us to use the same shock-capturing code to investigate the 
dependence on the scale size of the problem without introducing ad-hoc numerical stabilization techniques. This is in 
clear contrast to the EMTF model where the displacement current is retained in the Maxwell equation. In the EMTF 
model, there remain high frequency waves (electromagnetic and Langmuir waves) even in the long wavelength limit. 
These waves impose a severe restriction on the time step of an explicit time integration scheme. We thus think, in 
situations where high frequency waves do not play a major role, our model is better than the EMTF model in practice. 

One may expect that the QNTF model gives a good approximation to the EMTF model when the following 
condition is met 


Or 


u>„ 


V/u 


<s 1, 


(28) 


where Q. ce = q e B/m e c is the electron cyclotron frequency and V a, e = B / -y/4 np e is the electron Alfven speed (~ 
sjmfm,. times the Alfven speed). In this case, the two time scales are well separated and the interaction between them 
may be assumed to be weak. The above condition thus implies that relatively slow (slower than the electron cyclotron 
period) time scale phenomena being described by the QNTF model are almost decoupled from higher frequency 
phenomena ~ a> pe . It is also possible to estimate a normalized charge density fluctuation amplitude using the Gauss 
law V ■ E = 4np c (where p c is the charge density) as: 



where L and V are typical spatial scale length and velocity, respectively. Note that the electric field strength is 
estimated by E ~ VB/c. If one takes L ~ c/oj pe and V ~ Va, c , the right-hand side becomes O 2 e luf pe . This again 
suggests that the QNTF model is appropriate when Eq. (l28l ) is satisfied. 

From this analysis, we confirm that it is reasonable to neglect the displacement current for modeling low frequency 
non-relativistic plasma phenomena (i.e., VA.e/c <sc 1). Strictly speaking, however, it does not prove the validity of our 
assumption, which must be tested ultimately by direct comparison between the QNTF and EMTF models. This will 
be addressed in a future publication. 

We note that, even if the condition Eq. (l28l > is satisfied, application of the present model to phenomena of scale size 
less than the electron inertial length should be done with care. This is because electron kinetic effect is not properly 
taken into account in a fluid model, which will, however, play a role at this scale unless the electron fluid is unusually 
cold. Nevertheless, the inclusion of finite electron inertia is of critical importance at least for numerical stability even 
in the absence of kinetic effect. 


2.3. Effect of Resistivity 

In analogy with MHD, it is easy to understand that the parameter t) appearing in the generalized Ohm’s law 
Eq. (IT2l) represents the resistivity in the usual sense. However, at scales comparable to or less than the electron inertial 
length, it does not lead to magnetic diffusion. 

To demonstrate this, consider for simplicity a resistive medium so that the left-hand side of Eq. balances with 
the resistive term. For a long wavelength mode kclw pe <s 1, the induction equation for the magnetic field becomes a 
diffusion equation 

d ric 2 1 

* B * ^ V ‘ B ’ (30) 

from which one can immediately understand that 77 is actually proportional to the magnetic diffusivity. This is consis¬ 
tent with the usual resistive MHD equations. For kc/a> pe » 1, on the other hand. Ohm’s law becomes 

rVxVxE^^AJ. (31) 
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In this case, we obtain the following equation by taking rotation of the induction equation 


a. A i 

Hi 1 ” 


( 32 ) 


This is a pure damping equation that does not involve any spatial derivatives. In contrast to the diffusion equation, the 
electric current is damped locally without propagating to neighboring regions in this parameter regime. Therefore, the 
restriction on time step required for numerical stability, which would be severe for a diffusion equation in a strongly 
resistive medium, is very much relaxed. This form of resistivity was recently used in Amano et ah (2014) to improve 
the stability of a hybrid code in and around a vacuum region. 

This difference may be understood as follows. Physically, a finite // arises if there is a friction between the two 
fluids. The friction may be either due to Coulomb collisions or wave-particle interactions associated with unresolved 
microscopic turbulence (i.e., anomalous resistivity). For scales less than the electron inertial length, neither electrons 
nor ions are frozen-in to the magnetic field line. Therefore, the collision does not directly alter the magnetic field 
evolution. On the other hand, the relative streaming between the two fluids decays exponentially due to the friction, 
so does the electric current. We should mention that although the resistivity works differently in the two different 
parameter regimes, in both cases, the system relaxes to the same current-free state V x B = 0. 

In the following, the normalized resistivity defined as 


* _ w p 

11 =4^ 


(33) 


is used. Here or p = A = £ s cj 2 ps is defined with the local density, and thus, the collisionality is assumed to depend on 
the local plasma frequency. This is reasonable for modeling anomalous resistivity as its effective collision frequency 
will be characterized more or less by the plasma frequency. 


3. Numerical Algorithm 


We now describe the numerical algorithm used in our newly developed 3D simulation code. As has already 
become clear in the previous section, the QNTF equations consist of two coupled subsystems: conservation laws 
for five scalar (conservative) variables Eq. Q and the induction equation for the magnetic field Eq. (0J. Numerical 
solutions must satisfy the divergence-free condition V • B = 0 as much as possible so as to minimize the loss of 
accuracy in multidimensions (e.g., Balsara & Spiceil 1999 : Tothl 200d) . 

Because of the formal similarity with the MHD equations, we can apply some of numerical methods devel¬ 
oped for MHD. Here we adopt the HLL approximate Riemann solver combined with the UCT scheme (HLL-UCT; 
Londrillo & Del Zanna. 200(l 2004 ). that satisfies the divergence-free condition up to machine accuracy. The scheme 


was originally developed for classical MHD, and successfully applied to relativistic MHD as well (Del Zanna et al 


2003, 2007). Although the generalized Ohm’s law used in this study differs considerably from the ideal MHD, we 
find that the concept of the UCT is still useful for the system considered in this paper. 

Below we discuss only approximation in space; i.e., the temporal derivative always remains analytic. In this paper, 
we employ a scheme with second-order spatial accuracy for which finite difference and finite volume discretizations 
are identical. However, since it is well known that finite difference is computationally more efficient in multidimen¬ 
sions, we describe the numerical method with finite difference representation. This makes it easy to extend the scheme 
to higher orders if desired. The semi-discrete form of equations may be integrated using any stable time integration 
schemes^ In the numerical examples shown in section 4, we always adopt the third order TVD Runge-Kutta method 
of Shu & Osher ( 1988 ). 


3.1. Discretization in Space 

Let us consider a cartesian uniform mesh with sizes Ax, Ay, A z in each direction. We define the fluid conservative 
variables U (as point-value representations) at cell centers. In contrast, three components of the magnetic field vector 
B x , B v , B- are defined at staggered locations. Namely, each component of the magnetic field is defined at face centers 
along the normal direction (say, x direction for B x ). Here by a face center we mean the center of the two-dimensional 
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(2D) plane that defines the interface between neighboring cells in a specific direction. The same staggered discretiza¬ 
tion of the magnetic field is often employed both in PIC and MHD simulation codes. The reason for this is that the 
divergence-free condition for the magnetic field can be automatically preserved up to machine accuracy when the 
induction equation is integrated in time using the numerical flux (electric field) defined at edge centers (the center of 
the line separating two neighboring faces). The staggering technique applied to MHD is specifically referred to as the 
CT method (E vans & Hawley. 1988). In the next subsection, we discuss how the magnetic field update accommodates 
consistently with an approximate Riemann solver used to advance conservative variables defined at cell centers. 


3.2. HLL-UCTScheme 

Let us first consider a one-dimensional (ID) hyperbolic conservation law: 

dt dx 


(34) 


Temporal evolution of the solution vector u, (i indicates the index for a cell) defined at cell centers may be described 
by the following equation in the semi-discrete form: 


d 1 /, 

a u ' + s; ('■ 


i+l/2 - f/-1/2 


2 ) = 0. 


(35) 


The numerical flux fi+ 1/2 is defined such that the above difference equation gives an approximation to the spatial 
derivative with the desired accuracy. Note that although this equation looks similar to finite volume discretization, 
they are generally different when higher than second order schemes are concerned. (In a finite volume scheme, u,- 
and f,+i /2 in the above equation should be replaced by the cell average u, and point value f,+ 1 / 2 , respectively.) This 
illustrates only a ID conservation law, but extension to multidimensions is straightforward as far as the semi-discrete 
form combined with finite difference discretization is employed. In contrast to this, the situation becomes much more 
complicated when finite volume discretization is used in multidimensions. In any case, since below we only consider 
a second order scheme, these two approaches are identical. Nevertheless, the difference must be kept in mind to be 
prepared for extension to higher orders. 

The key question is how to evaluate the numerical flux ff+i/ 2 - A typical strategy is to reconstruct the left and right 
states of the solution vector at the cell interface as point-value representations, which we denote by u R |;2 and u ; K _ | /2 , 
respectively. One may then solve the Riemann problem at the cell interface either exactly or approximately to obtain 
the numerical flux. In general, the exact Riemann solver is very expensive and usually approximate Riemann solvers 
are a dopted. Although a lot of Riemann solvers have been proposed o ver the decades for the ideal MHD equations 
(e.g., Brio & Wu . 19881: Balsara . 1998 : Mivoshi & Kusano . 120051 2008 ), most of them cannot be applied to the QNTF 
model as they rely on the eigenstructure of the basic equations. Furthermore, since the presence of dispersive waves 
having characteristic temporal and spatial scales is an intrinsic nature of the QNTF model, the solution to the Riemann 
problem is no longer self-similar, making the situation much more complex. 

To avoid complication inherent in the physical model, we adopt the HLL approximate Riemann solver which does 
not require eigenmode decomposition. The numerical flux in this approximation (as a point-value representation) is 
given by 


f = 


a + f L + a f R - a + a (u R - u L ) 


a + + a 


(36) 


where f L ' R = f(u L ’ R ) and o: ± represent the maximum characteristic speeds (defined as absolute values) in the positive 
and negative directions, respectively. In general, f and f are different and the correction must be taken into account 
for higher than second order schemes (e.g., Del Zanna et al. . 2003 , 20071) . However, in a second order scheme used 
in this paper, taking f » f is sufficient. Note that we omit indices for a 1 , but they must always be evaluated at the 
boundary where the Riemann problem is defined. 

The HLL flux is obtained by assuming that the physical state is constant over the Riemann fan. Thus, the only 
spectral information required for the solver is the expansion velocities of the Riemann fan, which are estimate d by the 


maximum characteristic speeds a*. The scheme is also known as the central-upwind scheme (Kurganov et al.. 20011) . 
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which is by construction free from characteristic decomposition. Note that when the symmetry over the Riemann fan 
a + = a~ is further assumed, it reduce s to the well-known LLF (local Lax-Friedrichs) flux. In this case, it is also 
referred to as the central scheme ( Kurganov & Tadmor . 2000 1. 

Now we consider application to the QNTF model with the CT-type discretization described in the previous sub¬ 
section. Although the CT scheme has been widely used with many different Riemann solvers as well as reconstruction 
techniques, it is not a trivial question how to couple the magnetic fields defined at face centers and cell centers. If one 
tries to apply a ID scheme to a multidimensional problem via dimension by dimension (either with dimensionally spit 
or unsplit) approach, the magnetic field must be defined at the same location as other conservative variables. One may 
then obtain numerical fluxes at face centers for each direction. On the other hand, the CT-type discretization requires 
the numerical flux (i.e., electric field) defined at edge centers. One immediately notices that since the electric fields at 
edge centers are not available from a ID Riemann solver used to solve the fluid conservative variables, interpolation 
is needed to obtain the electric field at edge centers. This leads to many different variants of met hods that have been 
proposed in the literature (e.g., Rvu et al. . 1998 : Dai & Woodward! 1998 : Balsara & Snicer . 1999). 

The UCT framework gives a consistent way to calculate the numerical flux at edge centers. It is actually designed 
to be consistent with an underlying Riemann solver used to compute numerical fluxes for fluid conservative variables. 
For simplicity, below we consider a 2D version of the scheme in the x—y plane, but extension to 3D is trivial. The 
induction equation may be written as 


d c / ~ ~ \ 

— B x -i+l/2,j - — [E z; i + i/2J+l/2 ~ E z - i+ i/2,j-l/2) = 0 , 


dt 


-B 

dt 


'y;i,j+1 /2' 


st( £ 


z;i+l/2,j+l/2 ~ E z -i-\/2J+ll2 


2 ) = 0 . 


(37) 

(38) 


This form clearly suggests that the numerical flux £*;;+ 1 / 2 , 7 + 1/2 must be defined in a genuinely multidimensional man¬ 
ner because it simultaneously provides the flux for B x in y direction and B y in x direction (see also, Gardiner & Stone! 
20051). 


Londrillo & Del Zanna ( 20041) proposed the following formula to calculate the electric field at edge centers 


E z = 


a x a yE z 


+ — T-lL'xRy 

■ a x a y E z 


— + j-,R X Ly 

a x a y E - 


_ _ pRxRy 

+ a x a y t z 


a. a,. 


By - B\ 


a y a y 


B* y 


B r 


(arj +a x )(a* 


a y ) 


at + a r 


a J + a y 


(39) 


where the superscript indicates left and right states in x and y directions, respectively. (For instance, L x represents 
the left state in x direction.) The maximum characteristic speeds in x and y directions respectively are denoted by 
a x and a y . Again, E z ~ E z is satisfied to second order, but correction must be taken into account for a higher 
order approximation. It is clear that the above flux formula involves four states rather than two in the ID HLL flux. 
Therefore, this provides a flux of fully 2D in nature. Indeed, this coincides with a 2D central-upwind (or HLL) flux 
formula for a Hamilton-Jacobi equation given by Kurganov et al. (2001). This is actually to be expected because the 
induction equation in 2D may be recognized as a Hamilton-Jacobi equation in terms of the vector potential. 

The crucial point in the above flux formula is that it automatically and correctly reduces to the ID HLL flux when 
homogeneity in one direction is assumed. It is thus called the HLL-UCT scheme (Londrillo & Del Zanna, 2004). The 
second and third terms play the role for the upwind property, but had been ignored in earlier attempts to combine 
Riemann solvers with the CT discretization considering only averaging of the electric field (i.e., the first term). Notice 
that more elaborated multidimensional Riemann solvers have been presented recently in the literature (e.g., Balsara, 
20ic| 2012 ). In this respect, the HLL-UCT may be regarded as one of (and the simplest version of) the 2D Riemann 
solvers. 


Although the original work seems to be using implementation specific to the ideal MHD Ohm’s law (see, Del Zanna et al 


2007), application of the UCT scheme should not be restricted to a specific form of Ohm’s law. Indeed, once the 
electric field is determined at cell center (see section [3~4l) . the numerical flux £.- ; ,+ i/2,/+ 1/2 may be obtained by an 
appropriate reconstruction. However, the above form is not necessarily convenient as it formally involves 2D inter¬ 
polation of the electric field to edge center points. In the next subsection, we introduce slightly different definition 
of the numerical flux without loosing its advantage, which can more easily be implemented using successive ID 
reconstructions. 
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3.3. Reconstruction 

One has to consider spatial reconstruction of physical quantities to calculate the numerical fluxes. Throughout 
this paper, we use a piecewise linear polynomial with the Monotonized Central (MC) slope limiter for non-oscillatory 
reconstruction, which is therefore second order in space. On the other hand, interpolation of the magnetic field defined 
at face centers (i.e., primary data) to cell centers is performed as follows 


‘x;i+l/2.j + B x 

;+1/2-r) 

(40) 

'y,i,j+ 1/2 + B y, 

i.j- 1/2) - 

(41) 


which is also correct to second order accuracy and does not involve a nonlinear slope limiter. We note that there must 
be a correction term when one employs the divergence-free reconstruction technique as proposed by Balsam ( 200li 
2004). Nevertheless, we have not implemented it yet at present. 

We take (p,, v,, p,-,p e , v e , p e , E, B) defined at cell centers as the variables for which the reconstruction is performed. 
Recall that the density and velocity of ions and electrons are not independent quantities. However, since the conversion 
from ion to electron velocity involves calculation of V x B (i.e., spatial derivatives), it is rather convenient to determine 
the electron velocity first at cell centers which is then mapped to face centers via reconstruction for calculation of the 
numerical flux. Although this potentially introduces an inconsistency between ion and electron quantities (which 
must be related by the constraints Eqs. l[6][7])) reconstructed at face centers, we have not encountered any difficulties 
associated with this. For the same reason, we perform reconstruction for the electric field defined at cell centers, for 
which one needs to solve the generalized Ohm’s law (see section [3~4l >. 

Now consider calculation of the numerical flux F A at each x-face for the fluid conservative variables U. This 
requires the left and right states of fluid primitive variables p,, v,-, p,,p e , v e , p e and the transverse components of elec¬ 
tromagnetic field E y , E~, By, B-. Note that one can use the normal component of the magnetic field B x already defined 
at this point without any reconstruction (hence no ambiguity), whereas the electric field E x is not needed in the flux 
calculation. One now obtains the numerical flux F v using the reconstructed left and right states and the HLL flux 
formula Eq. (l36l >. The same procedure is applied in the y direction to obtain F v . 

In computing the numerical flux F^, one also calculates E. at each face as an appropriate average using available 
reconstructed left and right states. Writing the reconstruction procedure to obtain the left and right states symbolically 


we use the following HLL average 

{E--jj)i +\/2 

(E Zy ij)j+\/2 


fL,R . 

J 1+1/2 - 



(42) 

a x^f+ 1/2 1 

l Ez;i,j) + °T^i+ 1/2 ( 

E z-,i,j) 

(43) 


a+ + a~ 


<^Il/2 

(E z ,j) + a ^l 1/2 


(44) 


(*y + a- 

’ 


for the averaged electric field defined at x and y faces, respectively. Here the angle bracket () indicates the HLL 
average along the direction specified by the subscript (i + 1/2, j +1/2 respectively indicate averaging over x and y 
directions.) 

These HLL-averaged electric fields are further reconstructed and averaged in the other direction. Consequently, 
the electric field defined at edge center E z -j+i/ 2 ,j+ 1/2 is obtained as 


E z ;i+l/2,j+l/2 — 2 [((Ez;i,j)i+l/2)j+i /2 + j+ 1 / 2 ) 1 + 1 / 2 } 

■ _ nR x _ ijLx 

a x a x y;i+l/2,j+l/2 D y\i+l/2,j+l/2 


OtyOy 


„Ry _ nL y 

D x-,i +1 /2,j+ 1/2 " x-i +1 /2,j +1 /2 


az + a r 


(45) 
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where 


B 


Ly.Ry 

x;i+\/2,j+l/2 


rt^Ly,Ry 

*j+ 1/2 


(B x -MU2,j) > 


r L x ,R x 

D y;i+l/2,j+l/2 



y,i,j+ 1/2 


) 


(46) 


are results of ID reconstruction of the magnetic field. Here the first term represents the arithmetic mean of successive 
ID reconstruction-averaging procedures. For instance, the first term in the curly bracket of Eq. (l45l) represents the 
averaging in the x direction followed by y direction. It is readily be seen that this numerical flux reduces to the original 
definition Eq. ( l39l ) for the first-order piecewise constant reconstruction with constant maximum characteristic speeds. 
For a smooth region without discontinuities, one may also assume that a direct 2D reconstmction and successive 
ID reconstructions to the edge center point will give the same result within the order of accuracy of reconstruction. 
Although they do not necessarily coincide in general when a nonlinear non-oscillatory reconstruction is used, it is 
rather important that the above numerical flux retains the upwind property. It is again easy to confirm that by taking 
the ID limit this definition also reduces to the correct ID HLL flux. Therefore, we think it takes into account the 
essential feature for the UCT scheme. Our numerical experiments performed with Eq. (l45l) actually support this 
argument. 


3.4. Calculation of Electric Field 

The generalized Ohm’s law used in this study is written in an implicit form. Therefore, when discretized on a 
mesh, one must solve a matrix equation to obtain the electric field. Notice that once the magnetic field and fluid 
moment quantities are given, the equation is linear and thus solvable using any matrix solvers. As we have already 
seen in section lT2l we can reasonably assume V • E « 0 for low frequency phenomena in a non-relativistic plasma. 
We thus adopt the approximation 


VxVxE* -V 2 E 


(47) 


throughout in this paper. 

We define all the three components of the electric field at cell centers. Since we employ a second order scheme for 
solving the basic equations, second order central difference discretization is sufficient for approximating the spatial 
derivative in the generalized Ohm’s law. Thus, the left-hand side of the equation may be written as 

(A£ - c 2 V 2 S'j, * (a ijjf + 2(e x + e y + e z )) Gij.k 

~ £x (Ci— I Gi+l,j,kj ~ £y (fiij-l ,k G/./+ | ^ j — 6- (&i,j,k~\ EiJ.k- (-1^ (48) 

where e x - c: 2 / Ax 2 , e y = c 2 /Ay 2 , e, = c 2 /Az 2 and £ = E x ,E y ,E- represents the three components of the electric 
field vector. Likewise, the second order central difference is used to evaluate spatial derivatives in the right-hand 
side, which gives the source vector to the matrix equation. The discretization of the source term may possibly cause 
even-odd decoupling, but so far we have not found any problems. 

The resulting matrix equation is iteratively solved using the symmetric Gauss-Seidel method. Although it is known 
that the convergence of the method is not so fast in general, it is sufficient for our purpose because the matrix equation 
is diagonally dominant for the parameter regime of our interest. More specifically, the ratio between the diagonal and 
non-diagonal coefficients is proportional to /i 2 /(c 2 /w 2 e ) where h represents the grid size. Therefore, as far as the scale 
size large compared to the electron inertial length is concerned, it is essentially a diagonal matrix and the inversion is 
trivial. This is natural because the generalized Ohm’s law is no longer implicit in the absence of finite electron inertia. 
In most of numerical examples shown in section 4, the grid size is comparable to or larger than the electron inertial 
length, for which the matrix is relatively easy to invert. The convergence is checked by monitoring the residual during 
the iteration. It converges within a few iterations for a normalized tolerance of 10 2 in most of the problems. 

In cases where the electron inertial length is resolved by several grid points, sometimes 30 or more iterations are 
needed. If one is interested in the electron scale physics, it is desirable to adopt a more advanced numerical method 
such as multigrid. Nevertheless, our primary motivation is to introduce a physical upper bound to the phase speed of 
whistler waves for better numerical stability. For this reason, we have not yet tried to implement such complicated 
methods. 
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3.5. Characteristic Speed 


The only remaining task is to evaluate the maximum characteristic speed required to determine the width of 
the Riemann fan. In principle, it may be obtained by solving a linearized eigenvalue problem at the cell interface. 
However, it will involve much more cumbersome algebra than for the ideal MHD case. Furthermore, non-self- 
similarity of the problem makes such an approach impractical. Therefore, we here adopt a simplified approach as 
explained below. 

We first evaluate the wave phase speed with respect to the plasma rest frame by performing standard linear analysis 
of the basic equations for a homogeneous plasma (see. Appendix B for detail). After tedious but straightforward 
algebra, one finds that a wave phase speed X normalized to the Alfven speed Va = B/ yjAnipi + p,.) may be obtained 
by solving the 6th order polynomial equation 


PX 6 -QX 4 +RX 2 -S = 0 


(49) 


with the coefficients defined as follows 

P = (l+ar)\ (50) 

Q=1+/3 + (1+k 2 ) cos 2 6 + etc (l - cos 2 6 + 2 /?) + e 2 /c 2 (cos 2 6 + /3k 2 ^ , (51) 

R = (l +(2 + k 2 )P + s 2 k 2 /3)cos 2 6, (52) 

S - f3 cos 4 6. (53) 


Here s = m e /m.i, [3 = Anyipi + p e )/B 2 and 8 represents the wave propagation angle with respect to the ambient 
magnetic field. Note that the wavenumber k is normalized to £i n / Va which differs from the inverse ion inertia length 
by a factor 1 / Vl + s. It is easy to confirm that the characteristic equation reduces to the ideal MHD dispersion 
relation by taking the long wavelength limit k — > 0. A finite wavenumber k 4- 0 thus describes the two-fluid effect. 
The dispersion relation of Hall-MHD is obtained in the limit s —> 0 (i.e., negligible electron inertia). 

In the following, we always evaluate the phase speed at the Nyquist wavenumber. This choice is motivated 
from the conjecture that the Riemann fan should be bounded by the maximum possible characteristic speeds, which 
typically correspond to the wave at the Nyquist (i.e., largest) wavenumber. We note this does not hold when the grid 
size is much smaller than the electron inertial length. Nevertheless, we have confirmed that a different choice of the 
wavenumber does not introduce appreciable difference in our numerical results. 

A root with the maximum absolute value of the above equation gives the maximum wave phase speed in the plasma 
rest frame. A reasonable estimate for the expansion speed of the Riemann fan may be obtained by transforming it 
back to the laboratory frame. For this, one has to find the “bulk velocity” with which the bulk fluid moves with respect 
to the laboratory frame. This is, however, not a trivial problem in the present model. Namely, in situations where a 
finite current flows, the ion and electron flow velocities are different. The effect is more and more pronounced where 
the ion and electron dynamics is decoupled and the Hall current plays the role. This ambiguity arises perhaps because 
our perturbation analysis is performed for a homogeneous plasma where the current is assumed to be zero. Again, our 
strategy is to estimate an upper limit of the characteristic speeds. For this, we adopt 


a* = max {XVa ± v,-, XVa ± v e , 0) 


(54) 


as the maximum characteristic speeds, where v, and v e are ion and electron bulk velocities normal to the interface. 

The characteristic equation is cubic with respect to X 2 , which can easily be solved by an iterative method. Note 
that since we can use the solution obtained at the previous time step as the initial guess (which is likely to be very 
close to the solution), iterative methods would be superior to analytical calculation. We use the standard Newton 
method for root finding. It is easy to confirm from elementary calculus that the desired root should be found in the 
range X 2 >(Q+ Ve 2 -3PR)/(3P). 

Our approach is obviously not exact but gives a measure for the upper bound to the expansion speed of the Riemann 
fan. This implies that numerical dissipation is slightly increased with this approach. Nevertheless, the ambiguity in 
the definition of the Riemann fan appears only when the dispersive effect becomes important. In this case, a shock is 
no longer a pure discontinuity but has a structure of finite width accompanied by dispersive waves. For this reason, 
we think that the increased numerical dissipation is not a critical issue in this model. 
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3.6. Summary of Numerical Procedure 

Here we summarize the procedure used in our simulation code. We initialize the fluid primitive variables and 
the magnetic field at edge centers, such that V • B = 0 is satisfied. The cell-center magnetic field is calculated by 
interpolation, and the primitive variables are converted to the conservative variables at cell centers. 

At the beginning of the Runge-Kutta integration, the maximum characteristic speed is calculated at cell centers, 
and is assumed to be constant during the one Runge-Kutta time step. Appropriate mapping of the characteristic speed 
from cell centers to face/edge centers is performed. We take the maximum of the two neighboring characteristic 
speeds in the normal direction to obtain the estimate at each face. Then, their simple arithmetic average in the 
transverse direction is taken to evaluate the values at edge center. 

With this preparation, the following six steps are performed for each substep of the Runge-Kutta time integration: 

1. The electric field is calculated at cell centers by solving the generalized Ohm’s law Eq. (El). 

2. Reconstruction of fluid primitive variables and transverse electromagnetic field components defined at cell 
center is performed to obtain the left and right states at face centers for each direction. 

3. Numerical fluxes for the fluid conservative variables are calculated using the HLL flux formula Eq. (l36l) at each 
face. At the same time, the HLL average of the transverse electric field components are calculated according to 
the definition Eas. (143H44l> and are stored on a working array. 

4. Reconstruction of the electromagnetic field (i.e, the primary magnetic field and HLL-averaged electric field) 
defined at face center is performed. The numerical fluxes for the induction equation are calculated using the 
formula Eq. (l45l i. 

5. The fluid conservative variables at cell center as well as magnetic field at face center are updated in a conserva¬ 
tive manner using the numerical fluxes obtained above. 

6. The updated magnetic field at face centers is interpolated to cell center. The fluid primitive variables are then 
recovered at cell center from the updated conservative variables. 


This completes the description of the numerical procedure. Since our numerical algorithm follows earlier studies using 
the HLL-UCT with finite difference framework ( Del Zanna et all 20031 2007t) except for implementation specific to 
MHD, it is relatively easy to extend the code to higher orders. Implementation of higher order schemes is left for 
future studies. 


4. Numerical Results 

In this section, we discuss numerical results for several benchmark problems obtained with our new simulation 
code for the QNTF equations. Since we think it may replace the role of Hall-MHD and/or EMTF equations in 
plasma simulations, it is fair to employ benchmark problems that have been well tested with these models. We have 
found that, however, test problems for Hall-MHD models have not necessarily been well established yet; most of 
previous studies tested their codes against simple whistler wave propagation and the well-known GEM (Geospace 
Environment Modeling) magnetic reconnection challenge problems. In particular, problems involving shocks and/or 
other discontinuities in multidimensions have not been tested in the published literature. This may be related to the 
fact that many Hall-MHD codes implement numerical dissipation for short wavelength modes in an ad-hoc manner 
for numerical stability. 

We thus adopt benchmark problems that have been tested with the EMTF and MHD models. Since the QNTF 
model correctly reproduces the MHD result if the grid size is appropriately chosen, benchmark problems for MHD can 
be naturally extended to the regime where two-fluid effect becomes apparent. Concerning numerical results obtained 
with the EMTF model, we expect the difference should be negligible whenever the charge neutrality assumption is 
adequate. Therefore, for low-frequency {a> a> pe ) problems, it is reasonable to compare our numerical solutions with 
the published EMTF results. 

Since we know the maximum characteristic speed, the time step is chosen such that it always satisfies the CFL 
(Courant-Friedrichs-Lewy) stability criterion in the entire system. In the numerical examples shown in this paper, we 
always use a time step with a CFL number less than 0.5. More specifically, when the CFL number becomes greater 
than 0.5 (smaller than 0.25), the time step is halved (doubled). 
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Figure 1: Total mass density for the Brio-Wu shock tube problem at t = 0.1. Simulation results with three different resolutions are shown with 
black {N = 400), blue {N = 1600), and red (N = 6400) lines, respectively. 


In the numerical examples shown below, we use a polytropic index of y — 5/3, resistivity q - 0, and temperature 
ratio t - 1 unless otherwise stated. Since by construction our scheme automatically satisfies the solenoidal condition 
for the magnetic field, the V • B errors are not shown explicitly. We have confirmed that the error is indeed always at 
a level of machine accuracy. Although our code is fully 3D, only ID and 2D simulation results are discussed below 
because they are sufficient to demonstrate the capability of the code with modest computational time. 


4.1. Brio-Wu Shock Tube 

The first test problem is an extension of ID Riemann problem proposed by Brio & Wu ( 1988i) that is one of the 
standard test problems for the ideal MHD . It has also been tested by the EMTF mode l, which clearly demonstrates 
the appearance of dispersive features (e.g., Hakim et all 2006 : Kumar & Mishra . 2012 ). 

At the initial condition, the simulation box of unit length 0 < x < 1 was divided into the left and right states at the 
center x = 0.5. Their physical quantities were specified respectively as follows 


' p ' 


f L0 ) 

p 


1.0 

B.x/ V4tt 


0.75 

.By/ V4 K, 

lpft 

. 1.0 


' P ' 


( 0.1 f 

P 


0.125 

BJ V47T 


0.75 

.By/ V47T, 

right 

.-1.0 J 


(55) 


where p = p, + p e and p = /?, + p e are the total mass density and total pressure, respectively. Other quantities were 
zero everywhere at the initial condition. This setup is identical to the original Brio-Wu shock tube problem for the 
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Figure 2: Closeup view of ion and electron velocities for the Brio-Wu shock tube problem at t = 0.02. The x (top) and y (middle) components of 
ion velocity and y component of electron velocity (bottom) are shown. The meaning of the color is the same as the previous figure. 


ideal MHD. However, the parameter qi/nij can be chosen arbitrarily to control the ion inertial length A, oc (q l /m j ) 
Clearly, the ideal MHD corresponds to the limit <?,/»;, —> oo (A, —> 0). Here we chose q,lm, such that A, = 10 3 / \[p. 
In this case, the problem is close to MHD but two-fluid nature will appear when sufficiently small grid size is used 
to resolve A,. For this problem, y - 2 was used instead of 5/3. The ion to electron mass ratio was chosen to be 
rn.itm e = 100 . 


Fig. Q] shows snapshots of the total mass density p at t = 0.1 obtained with three different resolutions N = 
400,1600,6400 (i.e.. A* = 2.5 x 10 _3 ,6.25 x 10~ 4 ,1.5625 x 10 4 ). At the lowest resolution (black), the grid size 
was comparable to or larger than the ion inertial length and dispersive nature was not clearly visible. Consequently, 
the numerical solution was almost the same as the MHD. In contrast, a clear dispersive wave train structure appeared 
behind the slow shock at around x ~ 0.65 in the medium resolution run (blue). Similar structures were also found 
at the leading edge of fast mode rarefaction propagating to the left (x ~ 0.3), as well as around a compound wave 
(x ~ 0.5). These features became further pronounced in the solution with the highest resolution. It is interesting that 
such dispersive waves are clearly visible even with such a small ion inertial length. We found that the propagation 
speed of the slow shock appear s to be slower than the MHD result. This seems to be consistent with the result 
presented by Hakim et~akl (2006). This may be understood to be due to a finite Poynting flux carried by the dispersive 
waves behind the shock. 

It is noted that at this time (f = 0.1) whistler waves emitted to the right had already reached to the right-hand 
boundary, which is, however, barely visible in the mass density profile. The propagation of whistler waves can more 
clearly be identified in the velocity profiles at earlier times. Closeup view of the ion and electron velocities at t — 0.02 


16 






















0.10 

0.05 

cq °- 00 

-0.05 

- 0.10 


0.10 
0.05 

cq 1 0.00 

-0.05 

- 0.10 

0 10 20 30 40 50 60 

x 



Figure 3: Profiles of the in-plane wave magnetic field component obtained by the circularly polarized wave problem. The highest frequency case 
kAj ~ 2.24 (top), and lowest frequency case kAj « 0.22 (bottom) are shown respectively. Results for four different resolutions are shown with red 
(N = 16), green (N = 32), blue (N = 64), and black (N = 128) lines, respectively. Note that the blue and black lines almost coincide in this figure. 
Likewise, the exact solutions are indistinguishable from the highest resolution results. 


shown in Fig. [2] displays the propagation of whistler waves in both positive and negative directions again for three 
different resolutions. Note that the normal component of ion and electron velocities are identical in ID because 
(V x B) x = 0. 

We see that as increasing the spatial resolution, the leading edge of whistler waves became more and more ex¬ 
tended. In particular for the whistlers propagating to the right, the wavelength of the transverse electron velocity 
oscillations decreased in the positive x direction. This is clearly due to dispersive nature of whistler waves. Namely, 
since the whistler waves have dispersion relation to oc k 2 for Q,.,| «: to <sc |Q„,|, the group velocity becomes higher 
for shorter wavelength modes. Note that the maximum wave propagation velocity (both group and phase velocity) 
was limited by finite electron inertia effect in the highest resolution run because the grid size was comparable to the 
electron inertial length for this simulation ( nij/m e = 100). If the spatial resolution is insufficient to resolve these 
small-scale features, numerical dissipation automatically damps out these waves. 

Here, we have actually confirmed that our numerical simulation code automatically reduces essentially to ideal 
MHD when the grid size is larger than the ion inertial length. Since high frequency waves such as electromagnetic and 
Langmuir waves are absent, numerical stability only requires the CFL condition with respect to MHD characteristic 
speeds. We think this is a crucial advantage over the EMTF model. 
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Figure 4: Convergence of LI error for the in-plane wave magnetic field component as a function of the number of grid N for the circularly polarized 
wave problem. Convergence for three different wavenumber cases is shown with different colors: red for kA : x 2.24, green for kA, = 1.12, and blue 
for kAj » 0.22. The black dashed line indicates the theoretical second order convergence for reference. 


4.2. Circularly Polarized Wave 

It has been well known that a finite amplitude circularly polarized Alfven wave propagating along the ambient 
magnetic field (Bo) is an exact solution to the ideal MHD equations (e.g., Goldstein, 1978). Therefore, it is commonly 
used as a benchmark problem to test the accuracy of a numerical scheme. 

Generalization of the MHD exact solution to the QNTF model is possible. Actually, the dispersion relation of a 
finite amplitude circularly polarized electromagnetic wave is identical to the corresponding linear dispersion relation 


/ kOpi CO / Op e CO 

\ kc I o + \ kc I o + Q ce 


= 0 , 


(56) 


where Q, v = q s Bo/m s c is the cyclotron frequency. The positive (negative) o in the above dispersion relation indicates 
right-handed (left-handed) polarization. Adopting a coordinate system with orthogonal unit vectors e, (i = 1,2,3) 
such that ei is parallel to the background field, ei is a vector perpendicular to ei, and e 3 = ei x ea, we can write the 
eigenvector corresponding to the solution of the dispersion relation as follows 


B 2 = i : B() cos (kx - cot), B ] = i'Bo sin(fcr - ot), 

v Si 2 = V s cos (kx - ot), v s> 3 = V s sin(fct - ot). 


(57) 

(58) 
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where 


Vs = -f- 


-Q 


Q rv A: 


(59) 


The amplitude of the wave ( normalized to the background magnetic field Bq is a free parameter that can be chosen 
arbitrarily. Nevertheless, since the finite amplitude wave may suffer a parametric instability (e.g., Goldsteinl 1978 : 


Terasawa et al., 1986), we chose a small value /■' = 0.1 to suppress possible growth of the instability during the 


simulation. 

We set up a 2D simulation box with the computational domain 0 < x < 2L, 0 < y < L with the periodic boundary 
condition in both directions. The simulation box was resolved by 2 N x N grid points (i.e., Ajc = Ay). The background 
magnetic field was inclined with respect to the x axis by an angle 6 = tan ] (2). The wavenumber of a finite amplitude 
circularly polarized wave was chosen as ( k x , k y ) = (n/L, In/L) (i.e., mode = 1 in each direction). The initial condition 
was set up by rotating the exact solution by an angle 6 such that ei is parallel to the background magnetic field and 62 
is in the simulation plane. Since the wave propagates oblique to the grid, this provides a benchmark problem to test 
the accuracy of handling smooth profiles in multidimensions. In this simulation, time and length were normalized to 
the ion cyclotron frequency Q r , and ion inertial length A, = c/a> p j , respectively. The ambient magnetic field strength 
was chosen to be Bq — Vijr. 

By appropriately choosing L/Aj, one may test the code performance for different parameter regimes. We here 
discuss three different physical wavenumber cases, kAj ~ 2.24,1.12,0.22 (corresponding to L/A, = n,2n, 107r, re¬ 
spectively). The wave frequencies were respectively determined as w/Q„• « 5.72,1.89,0.25 by solving the dispersion 
relation for right-handed polarization (o> > 0). The lowest frequency case is essentially an Alfven wave in the MHD 
regime, whereas higher frequency cases correspond to dispersive whistler waves. In all cases, a plasma beta of /? = 0.1, 
and an ion to electron mass ratio of nij/m e = 256 were used. 

Rg. m shows the profiles of the in-plane component of wave magnetic field Bi = (-2 B x + B y )/ V5 as functions 
of x (i.e., cut through y = 0) after five wave periods. Results with the highest and lowest frequency runs are shown 
in top and bottom panels, respectively. The amplitude is normalized to the background field Bq. The simulations 
were performed with four different resolutions N = 16,32,64,128 for each case. We see that in both cases the 
numerical solutions converged to the analytic ones as increasing the resolution. Note that the initial conditions (or 
exact solutions) are indistinguishable from the highest resolution results and are not shown. 

The error convergence is shown in Fig. [4] for three different wavenumbers. The LI errors shown in this plot are 
those of the B 2 component evaluated as 


l UB 2 ) = 2 K,j(t = 5 T) - = 0)|, 


(60) 


where T denotes a wave period. This confirms that the code achieved roughly second order convergence (which is 
shown by the dashed line for reference). However, the error behavior became irregular for the higher frequency cases 
at high resolution. We have performed simulations with other wavenumbers and/or mass ratios, and confirmed that 
this irregular convergence appears when the grid size becomes comparable to or smaller than the electron inertial 
length. This may thus partly be related to the simplified estimate of the maximum characteristic speeds that is not 
strictly correct at this scale. In any case, since our primary focus is not on this small scale, we have not attempted to 
resolve the issue. 


4.3. Orszag-Tang Vortex 

The Orszag-Tang vortex problem has been one of the most standard benchmark problems for multidimensional 
MHD codes. It starts with a smooth initial profile, but the solution soon becomes complex involving many discon¬ 
tinuities. Therefore, it has been used to test the capability of handling multidimensional discontinuities which may 
introduce non-negligible magnetic monopoles in the numerical solution. On the other hand, to the authors knowl¬ 
edge, numerical results obtained with the Hall-MHD and/or EMTF models have not been published in the litera 


ture. Although there were numerical studies with similar but different initial conditions (e.g., Matthews et al., 1995 


Perri et al., 2012), they are not useful for comparison with the present results. Therefore, we generalize the original 
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problem to the parameter regime where the Hall term starts to play a role, but without changing the initial condition 
so that comparison with published results becomes easier. 

The computational domain was a unit square 0 < x < 1, 0 < y < 1 with the periodic boundary condition in both 
directions. The initial condition was given by the smooth profile 


v x = - sin (2ny), v y = sin ( 2nx ), 

B x = - Bq sin (2 ny ), B v = Bq sin (4nx ), 


(61) 

(62) 


where Bo = V47r. The density and pressure were constant; p = y 2 , p = y, and out-of-plane velocity and magnetic 
field components were zero v z — B z = 0. This definition is identical to the standard MHD test problem (except for 
normalization). Again, one needs to specify the ion inertial length (or qi/m ,) relative to the box size. 

Numerical results with four different initial ion inertial lengths A, = 10 3 , 10 2 ,5 x 10 2 , 10 1 are shown in Figs.0 
and [6] We used 200 x 200 grid points and m,/w e = 100 in all the simulation runs. Fig.[5]displays snapshots of the 
temperature distribution at t = 0.5. (Here, the temperature is defined with the total pressure and mass density as p/p.) 
In Fig. [6j the color indicates the out-of-plane component of the magnetic field normalized to Bq, whereas the black 
solid lines are contours of the vector potential A z (i.e., magnetic field lines in the plane) at the same instant of time. 
We confirmed that the V • B error in the numerical solutions was on the order of round-off error, which is consistent 
with the design of the numerical scheme. 

Since the grid size was Ax = 5 x 10~ 3 , the run with the smallest ion inertial length A, = 10 1 did not resolve 
the ion scale. It was thus essentially in the MHD regime and the results may be compared with published results 
for ideal MHD. We see that our simulation well captured discontinuities and other small scale features. This is 
not surprising because in this case our scheme becomes almost the same as the second order scheme described in 
Londrillo & Del Zannai (2004) except for subtleties. Although the overall structures were similar in all the cases, 
dispersive features clearly appeared in the numerical solutions as increasing the ion inertial length. The fact that 
the Hall term played the role in the numerical solutions can most clearly be visible in Fig. [6] because the out-of- 
plane magnetic field disappears in the ideal MHD limit. Even in the case with A, = 10~ 2 , in which the temperature 
distribution looks almost unchanged from the MHD result, amplitude of B z became substantial (at around x ~ 1.5,y ~ 
0.9 and its symmetric point). For larger ion inertial lengths, the solution became much more complex even in this 
early stage of evolution. 

The numerical results clearly indicate that it is possible to implement the two-fluid effect in a shock-capturing 
code without loosing its advantages. A shock-capturing code gives a non-oscillatory solution when encountered by 
discontinuities by automatically introducing required dissipation. On the other hand, dispersive character must be 
retained at ion and electron inertial scales, which tends to produce an oscillatory solution. Our numerical simulation 
code successfully implements these two contradictory features at the same time without suffering from numerical 
instability. 


4.4. Magnetic Reconnection 

Our final test problem is the GEM magnet ic reconnection challenge problem that has been well tested with many 


different simulation models (Birn et al., 2001). Here we use this problem to test the effect of re sistiv ity. The initial 
condition was given by the Harris equilibrium which is identical to the one described in Bimetal. (2001). The 
magnetic field and number density profiles were given by 


and 


B x (y) = Bo tanh (y/d) 


(63) 


n(y) = no sech 2 (y/d ) + nt, g , 


(64) 


respectively. Here, the width of the current sheet is represented by d. The ion and electron temperatures were 
determined by the pressure balance condition no(Tj + T e ) = B^/Sn. We here used a temperature ratio of r = 5 to be 
consistent with the original problem. The ion and electron drift velocities in the z direction must satisfy the relation 
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r = -ViyJv e j for the initial condition to be a Vlasov-Maxwell equilibrium. The drift velocities were thus initialized as 


cfio t no sech 2 (y/d) 

4ne 1 + r no sech 2 (y/d) + n bg ’ 
cB 0 1 «o sech 2 (y/flT) 

47 re 1 + tji 0 sech 2 (y/d ) + n bg 


(65) 

( 66 ) 


for ions and electrons, respectively. Other quantities were initialized with zero. 

The normalization of time and space was such that the ion cyclotron frequency O r , = eBo/m/c = 1 and the 
ion inertial length c/u) p i = c fnfAmioe 2 = 1. Accordingly, the velocity was normalized to the ion Alfven speed 
V A j = Bq/ '/‘XmiQin,. Note that in the normalized unit, Bq/ v/4/r = no = m, = 1. We used a rectangular simulation 
domain -L < x < +L and -L/2 < y < +L/2 with L = 12.8. The current sheet thickness was taken to be d - 0.5. The 
simulations were performed with 256 x 128 grid points with a mass ratio of mj/m e = 25. The grid size Ax = Ay = 0.1 
was thus roughly comparable to the electron inertial length for the initial current sheet density. The periodic boundary 
condition was used in x direction, while the conducting wall boundary was used in y direction. 

To initiate magnetic reconnection, the initial magnetic field was perturbed with the out-of-plane component of the 
vector potential given by 


where <J>o = 0.1. With this initial condition, the X-point was located initially at the origin from which the reconnection 
process started to evolve. 

Fig.|7]shows time evolution of reconnected magnetic flux for simulations with five different normalized resistivi¬ 
ties: rf = 0, 10 6 , 5 x 10 6 , 10 5 , 10 4 . Recall that since the resistivity is normalized with respect to the local plasma 
frequency, the diffusivity depends on the local density even with a spatially constant resistivity. The reconnected flux 
was computed by 


m 


i r +L 

— \B y (x,y = 0, t)\dx. 
2B o J-L 


( 68 ) 


In all cases except for rf = 10 4 , the reconnection rates estimated from the slope of reconnected flux exceeded ~ 0.1 
in the nonlinear phase, consistent with published results. 

Fig.[8]shows snapshots at Q n t = 21.5 for rf — 0, at which the normalized reconnected flux reached t/A 1 ) — 1.0. 
Fast outflows originating from the X-point are clearly seen for both ion and electron velocities. The ions and electrons 
were accelerated to the Alfven speeds defined for each species (V A j for ions and V Ae = fm/Jm^V A .i for electrons). The 
out-of-plane magnetic field B z was generated in association with the decoupling between ion and electron dynamics. 
Similar characteristics were also found in other runs at the time such that i//(t) = 1.0, although the magnitude decreased 
with increasing the resistivity. These features are consistent with previous numerical studies. 

There was a quantitative difference in the time evolution of reconnected flux in between f = 10 6 and 5 x 10 6 
after Q,,f > 25. We found that the further increase in the reconnection rate in smaller resistivity runs was associated 
with the formation of secondary magnetic islands. Figs. l9l and ITOlcomt) are runs with if - 0 and 5 x 10 6 at the time 
when the same amount of flux was reconnected i//(t) = 2.0 (Q„-f = 25.0 and 31.5 respectively). Two newly formed 
X-points in the elongated thin current sheet can be seen in Fig. [9] whereas the dynamics was still dominated by a 
single X-point in Fig. [TO] Formation of secondary magnetic islands was not observed until the end of simulations 
Cl c jt = 50 for rf > 5 x 10 () . This difference may be understood from the viewpoint of the stability of the thin (i.e., 
electron-scale) current sheet generated in the nonlinear phase of reconnection process. In higher resistivity runs, the 
length of the thin current sheet characterized by the fast electron outflow became shorter because of the diffusion 
during the propagation. This is probably the reason for prohibiting the excitation of a secondary tearing mode. The 


electron-scale current sheet may also become unstable against another instability in fully 3D simulations (Drake et al 


19941) . Nevertheless, one must bare in mind that electron kinetic physics, which is not included in the present model, 
will play a role in the dynamics of such a thin current sheet. 
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5. Conclusions 


In the present paper, we have proposed the QNTF model for collisionless plasmas. The basic equations are the 
fluid equations for ions and electrons and the Maxwell equations without displacement current. The absence of the 
displacement current indicates the charge neutrality, which is a reasonable assumption for low frequency phenomena. 
Rewriting the basic equations, it is shown that the system consists of the conservation laws for five scalar variables 
(total mass, momentum, and energy), and the induction equation for the magnetic field, supplemented by the general¬ 
ized Ohm’s law. The system fully takes into account finite electron inertia effect, which introduces the upper bound 
to the phase speed of whistler waves. On the other hand, it reduces to the ideal MHD in the long wavelength limit. 
No high frequency waves (such as Langmuir or electromagnetic waves) exist in the system, which would otherwise 
impose a severe restriction on the simulation time step even in the long wavelength limit. 

We have developed a 3D numerical simulation code for solving the QNTF equations. The code employs the 
HLL approximate Riemann solver combined with the UCT method for the induction equation. With this approach, 
we avoid complicated characteristic decomposition with keeping the advantage of a shock-capturing scheme. The 
UCT scheme guarantees the divergence-free condition for the magnetic field up to machine accuracy, without loosing 
upwind property. Therefore, the code is able to capture complicated multidimensional sharp discontinuities without 
appreciable numerical oscillation. At the same time, it successfully describes dispersive waves arising from the two- 
fluid effect without numerical problem. The upper bound of whistler wave phase speed appropriately introduced by 
finite electron inertia helps to overcome the well-known numerical stability issue in dealing with the short wavelength 
mode. 

In the numerical examples, we have confirmed that the numerical solutions reduce to the ideal MHD results when 
the ion inertial length is small enough. In this case, numerical stability requires only the CFL condition with respect to 
the MHD characteristic speed because of the absence of high frequency waves. This is in clear contrast to the EMTF 
equations in which high frequency waves must be resolved by the simulation time step. Since these waves will not 
play a major role in non-relativistic plasmas, the low frequency approximation is adequate. Consequently, we believe 
that the QNTF model offers a better alternative to the Hall-MHD and/or EMTF models. 
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Appendix A. Derivation of Generalized Ohm’s Law 


We here introduce a collision term in the right-hand side of the equation of motion for the two fluids to take 
into account finite resistivity rj. It is intended to be rather phenomenological (or anomalous) such as due to kinetic 
wave-particle interactions. The collision term is defined as 


R, - —R, 


77 le + U 2 pi 

4tt Qi 


qt q e 

—Pi\i + - p e X t 

rtii m P 


. ,2 

n T 

An J 


nti m ( 


mi m l 


(A.l) 


where R, and R e are for ion and electron fluids respectively. By adding it into the equation of motion, one may take 
into account effective friction between the species. It is easy to understand from the symmetry that the definition does 
not violate the momentum conservation law. 

The generalized Ohm’s law Eq. (Il2l) used in this paper may be obtained by taking a weighted sum between the 
two equations with the weight factor q s /m s . This yields 


| J + V ' 


Y — oav.s-v s + Ps i) = y 


q] F q 2 s v, 

Ps —E +Ps — — XB 


n 2 T 

- — OJ j. 
An p 


(A.2) 
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Using Faraday’s law, the first term may be rewritten as follows 



= V x V x E. 
4:t 


(A3) 


Now rearranging the equation, we arrive at 



which is identical to Eq. (IT2l) . It should be noted that contributions from both ion and electron fluids are fully taken 


into account in this form of Ohm’s law, so that it always provides the correct equation regardless of the ion-to-electron 
mass ratio. 

Appendix B. Linear Dispersion Analysis 


Here we outline the linear dispersion analysis for the QNTF model. We consider a homogeneous electron-proton 


plasma without resistivity, i.e., q, = -q,, = e and q = 0. The ion fluid equations, the generalized Ohm’s law, and 
the induction equation are used. Hereafter, we drop the subscript for the particle species unless necessary, and the 
fluid quantities without subscript must be read as the ion quantities. In addition, we use the definition for the thermal 
velocity Vf = yki-iTj/ni;, V 2 = ykBT e /m e (where is the Boltzmann constant) and the following notation: s = m e /mj, 
t = Ti/T e , /j — si 1 - ts)/(1 + s). Note that in the derivation below, the assumption of constant r is not necessary. 

It is easy to obtain the following equation from the linearized ion fluid equations 



(B.l) 


whereas for the generalized Ohm’s law, we have 



(B.2) 


with 



(B.3) 


(B.4) 


Rearranging the above equation using <5B = ck/w x SE, we obtain 



(B.5) 


nii \ c 7 nii 

where we have defined A 2 = V 2 /£lk = c 2 /u) 2 pi ( 1 + s) and £l cl = Q C! Bo/Bo, respectively. 

Now our task is to eliminate bv from Eqs. (IB. 1 1 ) and (IB. 51) . to get a dispersion matrix D satisfying 


D -SE = 0, 


(B.6) 


from which the dispersion relation is naturally obtained as |D| = 0. For this purpose, we introduce the following 
matrix notation: 


M = w 2 I - V^kk 

k 2 N • <5E = m ci x k x k x <5E - aojk x k x ()E 
W ■ 6\ = £l ci x 6\, 
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(B.8) 

(B.9) 










where the sound speed is defined by 


, , , yk K (T, + T e ) 

Vj = Vr + /uV 2 = r B e) 


nti + m e 


Solving the equations with respect to the ion velocity, we have 


8\ = -i—k 2 A 2 Mr x • N • 8E. 

Mj 

Substituting this into the generalized Ohm’s law Eq. (1B.5L we finally obtain 

V 2 1 

D = I + k 2 A 2 m ■ M _l • N + !i— kk • M _1 • N H— N 


(O 


U) 


as the dispersion matrix. 

Explicit form of the normalized dispersion matrix is given as follows: 


D = 


1 


X 2 (X 2 -P) 


[D mhd + D, + D e ], 


where 


D 


MHD 


X 2 - /3 2 )(X 2 - cos 2 9) 0 (X 2 - p 2 ) cos 9 sin 9 

0 X 4 - (P 1 + \)X 2 + p 2 cos 2 9 0 


0 


0 


D, = ikX 


I), = wkX(X 2 -p-) 


X 2 (X 2 -p 2 ) 

0 X 2 - p 2 + /up 2 sin 2 9 0 ' 

(X 2 -p 2 )cos 2 d 0 (X 2 -p 2 ) cos 6 sin 9 

0 /ip 2 e cos 9 sin 6 0 

-iKX cos 2 9 -1 iKX cos 9 sin 9 

cos 2 9 -ikX - cos 9 sin 9 

iKX cos 9 sin 9 0 -ii<X sin 2 9 


(B.10) 

CB. 11) 

CB. 12) 

(B.13) 

(B.14) 

(B.15) 

(B.16) 


In the above expression, X = a>/kVA is the normalized phase speed, k = kA is the normalized wavenumber, 9 is the 
wave propagation angle with respect to the ambient magnetic field. We also defined p = E 2 / V 2 and p e — V 2 /V 2 v 
respectively. 

One may easily understand that D mhd (which remains finite for k —» 0) corresponds to the MHD limit. On 
the other hand, D, (independent of s) and D e represent respectively the ion and electron inertia effects. Taking the 
determinant, and arranging it into the polynomial form, we obtain the dispersion relation given in Eq. ( l49l >. Most 
of cumbersome calculation presented in the derivation has been performed with a computer algebra system package 
SymPy ( SvmPv Development Teatrii 2014 ). 

Since the highest phase speed appears at the parallel propagation, we here investigate this special case in detail. 
For the parallel propagation, the sound mode decouples from the electromagnetic modes, and the dispersion relation 
may be factorized as follows 


(X 2 -P) {(1 + ek 2 ) 2 X 4 - (2 + (1 + £ 2 )i?)X 2 + l} = 0. 


(B.17) 


The dispersion relation for electromagnetic waves corresponding to the second factor (i.e., the curly bracket) is shown 
in Fig. [TT] It is clearly seen that the wave frequency approaches to the electron cyclotron frequency at the short 
wavelength limit. The phase speed, on the other hand, has a maximum at around k - yfmJnT e (i.e., kc/oj pi , — 1) and 
the maximum phase speed is given approximately by X ^ \jm ; /m,,/2. 

Note that, in a previous publication, we have shown that the deviation from the Hall-MHD dispersion a> oc k 2 
occurs approximately at k ^ (m,/m e ) l/4 (Amanpet ah, 2014, Appendix, B.). The weak dependence of the critical 


wavenumber on the mass ratio indicates that, at scale length on the order of or larger than the ion inertial length, the 
result should not depend strongly on the reduced mass ratio. 
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Figure 5: Temperature distribution for the Orszag-Tang vortex problem at t = 0.5. The solutions with different initial ion inertial lengths are shown; 
Aj = 10 -3 (top left), 10 -2 (top right), 5 x 10 -2 (bottom left), 10 -1 bottom right. 
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Figure 6: Out-of-plane magnetic field component B- and the vector potential A z for the Orszag-Tang vortex problem at t = 0.5. The color indicates 
B z , while black lines are contours of A z (magnetic field lines in the .v—v plane), respectively. The format is the same as the previous figure. 
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Figure 7: Time evolution of reconnected magnetic flux for the magnetic reconnection problem. Five runs with different normalized resistivities are 
shown in different colors. The black, red, green, blue, magenta lines correspond to rj* = 0,10 -6 ,5 x 10 -6 ,10 -5 ,10 -4 , respectively. 
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Figure 8 : Snapshots for the magnetic reconnection problem at Q cl -f = 21.5 for the run with rj* = 0. Shown are the mass density (top left), out-of- 
plane magnetic field component (top right), x component of electron velocity (bottom left), and ion velocity (bottom right). The white lines in each 
panel shows the magnetic field lines (contour lines for the vector potential A z ). At this time, the reconnected magnetic flux reached = 1.0. 



Figure 9: Snapshots for the magnetic reconnection problem at f l C jt = 25.0 (corresponds to the time = 2.0) for the run with 77 * = 0. The format 
is the same as the previous figure. 
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Figure 10: Snapshots for the magnetic reconnection problem at Q c ,-f = 31.5 (corresponds to the time = 2.0) for the run with rj* = 5 x 10 
The format is the same as the previous figure. 
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Figure 11: Linear dispersion relation for electromagnetic waves propagating parallel to the ambient magnetic field. Wave frequencies as functions 
of wavenumber are shown in the top panel for four different mass ratios rrii/m e = 25,100,400,1836. Corresponding phase speeds are shown in the 
bottom panel. Only the waves with maximum phase speed (i.e., on the whistler branch) are shown. 
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